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Hypersonic flow simulations using the node based, unstructured grid code FUN3D are 
presented. Applications include simple (cylinder) and complex (towed ballute) configura- 
tions. Emphasis throughout is on computation of stagnation region heating in hypersonic 
flow on tetrahedral grids. Hypersonic flow over a cylinder provides a simple test prob- 
lem for exposing any flaws in a simulation algorithm with regard to its ability to compute 
accurate heating on such grids. Such flaws predominantly derive from the quality of the 
captured shock. The importance of pure tetrahedral formulations are discussed. Algo- 
rithm adjustments for the baseline Roe / Symmetric, Total- Variation-Diminishing (STVD) 
formulation to deal with simulation accuracy are presented. Formulations of surface nor- 
mal gradients to compute heating and diffusion to the surface as needed for a radiative 
equilibrium wall boundary condition and finite catalytic wall boundary in the node-based 
unstructured environment are developed. A satisfactory resolution of the heating problem 
on tetrahedral grids is not realized here; however, a definition of a test problem, and dis- 
cussion of observed algorithm behaviors to date are presented in order to promote further 
research on this important problem. 


I. Introduction 

The vision for computational aerot her mo dynamic design and multi-physics analyses of complex configu- 
rations in hypersonic flow require simulation capabilities that are accurate (with quantifiable uncertainties) 
and robust. The descriptor robust has been used for many years to indicate subjective qualities of relia- 
bility; when a simulation is started it will always converge to a correct result if suitable initial conditions 
and boundary conditions are applied. For the purposes of this discussion, a robust scheme must also be 
adaptive - it must be able to start from any reasonable initial condition for surface and volume grids and 
automatically adapt the grid and algorithm until a simulation is obtained to user specified uncertainties in 
one or more design variables. A reasonable initial condition for grid is another subjective term implying that 
it may be quickly and automatically generated without a priori knowledge of flow topology or structural 
response to aerot her mal loads. 

Regardless of the underlying grid type (structured, unstructured mixed elements, unstructured Cartesian, 
overset structured or unstructured) a robust scheme must either: (1) be free from any requirement of special 
topological grid constraints in resolving all flow structures; or (2) be able to automatically adapt specially 
required grid topologies wherever they are needed in the flow field. In the context of simulation of hypersonic, 
blunt body, stagnation region heating it is noted that three-dimensional simulations using high aspect ratio 
tetrahedra to resolve the boundary layer and the bow shock produce poor heating. This statement is based on 
limited published discussion of this issue 1 4 and an inability to identify any other published source including 
a detailed examination of hypersonic heating on anything other than a semi-structured (prismatic) grid 
system in the stagnation region. 

Problems with heating can be overcome if one uses semi-structured grid (prisms) across the boundary 
layer and adapts the grid to the shock. Nompelis et. al. 3 show excellent heating results for the cylinder test 
problem 1 * on families of grids where a prismatic grid was adapted to the shock and acceptable heating results 

* Senior Research Engineer, Aerothermodynamics Branch; AIAA Fellow 
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when an unbiased (random face orientation), tetrahedral grid was adapted to the shock. Prismatic elements 
are relatively easy to generate on blunt body geometries and their superior performance with respect to 
aeroheating predictions have led to their use as standard practice in unstructured simulations. There is no 
better alternative with todays algorithms than use of prismatic elements to capture the boundary layer and 
bow shock. 

If one is doing a simulation without any free shear layers or internal shocks the specialized application 
of prismatic elements at the body and (possibly) bow shock is perfectly acceptable. If such flow structures 
exist - as they do in almost all interesting problems - then the accuracy of any algorithm that ignores special 
topological grid requirements where the features are harder to enforce must be questioned. This issue is not 
unique to unstructured simulations (though the consequences to heating are more severe in the unstructured 
world). NASA’s main computational aerothermodynamic simulation codes LAURA 5,6 and DPLR 7 both take 
great pains to align structured grid with the captured bow shock to improve solution quality but accept the 
limitation that internal shocks (e.g. wake recompression, control surface compression) are computed with no 
special grid consideration. In fact, sharp double cone code validation test problems 8 revealed that extremely 
fine structured grids were required to achieve a grid converged solution in a problem where structured grid 
could not be easily aligned with internally reflected shocks and shear layer over a separation bubble. 9, 10 

The position advocated herein is that unstructured grids provide the greatest flexibility to adapt to 
evolving flow structures (viscous and inviscid) and complex, deforming bodies in a hypersonic flow simulation 
without a requirement for significant user intervention. More simply put, unstructured elements provide the 
greatest opportunity to create a robust aerothermodynamic simulation. A discussion of the simulation needs 
for hypersonic flow over an inflated aerobrake in the next section will provide context for this position. 
The problem with the simulation of heating on tetrahedral grids is currently one of the biggest obstacles to 
realization of such a simulation capability. Two basic approaches to address this problem are identified. 

The first approach is to define an algorithm which produces consistently accurate results regardless of 
the grid element type. It may be speculated that quasi one-dimensional flux reconstruction is not up to the 
task here based on effort expended to date on this path. Finite element approaches may have something 
new to offer. 4 Grid adaption is likely needed as well but the additional requirement to form special elements 
aligned with the high gradient region is not considered part of this first approach. The second approach is to 
define an algorithm in which any specially required, high aspect ratio, prismatic elements are automatically 
generated where needed and they adapt and align as the shear layers and internal shocks evolve. 

The goal of this paper is to focus on the first approach to conquering this major obstacle - define an 
algorithm which produces consistently accurate results regardless of the topology of the cells. A robust 
solution to the problem has not yet been identified but the definition of a sensitive test problem and a review 
of algorithm options already tested are presented as an initial step to this end. 

II. Spacecraft - Tether - Ballute System Simulation 

Large, inflatable ballutes (balloon-parachutes) have been proposed as hypersonic decelerators for plane- 
tary aerocapture applications. 11 13 A simulation of a spacecraft - tether - ballute system using the unstruc- 
tured grid, Navier-Stokes flow solver FUN3D is presented in Figs. 1-2. These results, repeated from an 
earlier work, 14 are presented here to provide context for a position advocating use of an adaptive, tetrahedral 
grid system for hypersonic simulations. 

The ballute has a 52 m ring diameter and 13 m cross-sectional diameter. Conditions of the simulation are 
for a Titan Organics Explorer with velocity equal to 8550 m/s and density equal to 1.9 x 10 -7 kg/m 3 . 11, 12 All 
surface temperatures are set to a constant value equal to 500 K. The gas model includes molecular nitrogen 
and atomic nitrogen in thermochemical nonequilibrium. The towing spacecraft is a Pathfinder shape - 70 
deg spherically capped cone with a 6 m base diameter. The simulation domain encompasses a 90 degree 
wedge about the system axis. The spacecraft looks sharp at the nose in the figures because of the view 
angle combined with the 90 degree domain cutting planes. The simulation assumes symmetric flow with four 
0.3 m diameter compressive tethers attached to the toroid so that the leading edge of the tether is tangent 
to the toroid outer surface. A compressive tether is a flexible cylinder which can be inflated to withstand 
compressive loads and is used to position the toroid in space prior to entry. Tensile load tethers would likely 
be encased within the compressive tether. At present, only a continuum simulation is enabled. Flow over the 
compressive tethers is deep in the transitional flow domain in which the validity of Navier-Stokes analyses is 
inaccurate. Ideally, a fully coupled continuum - rarefied analysis would be brought to bear on this complex 
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system of disparate length scales. 

All of the flow features illustrated in Figs. 1-2 are simulated using only tetrahedral cells. The flow 
conditions are at a relatively low Reynolds number; consequently the aspect ratio of cells in the boundary 
layer does not need to be very large. There is a spacecraft bow shock that interacts with the toroid bow 
shock. The toroid bow shock creates a focussed compression in the center of the toroid that drives a reverse 
flow upstream toward the base of the towing spacecraft. A shock forms further upstream, diverting flow 
around the core. The compressive tether passes through the spacecraft bow shock and reenters the toroid 
bow shock. A merged shock layer about the tether itself interacts with these external features. The shocked 
flow off the tether spills out onto the ballute surface and a local hot spot is formed near the impingement 
point as seen in Fig. 2. The shear layer around the recirculating core regulates energy exchange between 
flow processed by the spacecraft bow shock and flow processed behind the spacecraft bow-shock and toroid 
bow-shock interaction. As Reynolds number increases with deeper descent into the planetary atmosphere the 
reverse flow moves further upstream and becomes unsteady. The simulation here assumes a rigid body but 
a robust simulation must also handle the multi-physics interactions associated with an inflatable structure 
that undergoes significant deformation. 

Obtaining a grid converged simulation on this problem is challenging. Obtaining the solution under 
constraints that prismatic grids must adapt and align to all of these developing features (shocks and shear 
layers) on deforming bodies is even more difficult. That is not to say such constraints could not be realized 
in a robust simulation capability - only that an algorithm without a requirement for such constraints is more 
likely to be robust. 



Figure 1. Unstructured grid used in simulation of spacecraft — tether — ballute system. Colors correspond 
to pressure levels nondimensionalized by pooV^. 


III. Numerical Tools 


A. LAURA 

Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) is a high fidelity, structured grid 
analysis tool, specialized for hypersonic re-entry physics, utilizing state-of-art algorithms for computational 
fluid dynamic (CFD) simulations. 5,6 Key elements of LAURA include Roes averaging 15 and Yee’s Symmetric 
Total Variation Diminishing (STVD) 16 formulation of second-order, inviscid flux. Yee’s STVD formulation 
has been found to be exceptionally robust and Courant-number-independent using first point-implicit and 
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Figure 2. Surface heating in W/cm 2 on ballute in vicinity of tether attachment point. 


then line-implicit relaxation for hypersonic flow simulations. The TVD algorithm uses a non-linear, minmod 
function as a flux limiter that maintains second-order accuracy away from extrema but can admit limit 
cycles in the convergence process, particularly in the vicinity of captured shocks. This occurrence usually 
manifests itself as a stalling of convergence at a very low error norm, essentially a benign ringing in the 
solution at a level that has no impact on aerothermodynamic quantities. Viscous flux is computed using 
central differences. 


B. FUN3D 

FUN3D is a node based, fully unstructured, finite volume solver of the Euler and Navier-Stokes equations. 17 
The method uses Least Squares (LS) gradient information to execute second-order accurate inviscid flux 
reconstruction. 18 Viscous gradients are computed from a Green-Gauss formulation. A suite of modules 
includes all of the gas physics models in LAURA and VULCAN 19 for thermodynamics, transport properties, 
chemical kinetics, and thermal relaxation. The inviscid flux reconstruction algorithms used in LAURA have 
been modified where necessary for an unstructured grid context and incorporated in the generic gas path 
through FUN3D. 1,2 A familiarity with formulation of conservation laws with Riemann solvers is assumed. 
Only elements of the algorithm that differ from the baseline FUN3D are presented here. 


1. Inviscid Flux Formulation 


The colored line segments in Fig. 3 indicate faces of the node based, median dual control volume used in 
FUN3D. The face that separates the control volumes about node L and node R crosses edge LR at midpoint 
M. The inviscid flux crossing face M for the generic gas path through FUN3D is defined in Eq. 1 using 
information available at the endpoints of edge LR. The flux £m is defined as the second-order accurate 
average of the left and right states of the edge minus a first- or second-order dissipation term that provides 
upwind bias to the reconstructed flux. 



|f L + i R - Rjv/IAmI [dq M - rfqjim]} 


( 1 ) 


The unit normal crossing cell face at point M is not generally aligned with the edge LR. The upwind 
dissipation is based on a quasi one-dimensional characteristic gradient. The characteristic direction is or- 
thogonal to the face at M. The eigenvector matrix R m that transforms from conservative to characteristic 
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Figure 3. Representation of stencil used for quasi one-dimensional reconstruction of inviscid flux. Dashed 
lines indicate dual volume. Solid lines indicate edges connecting nodes. 


variables is based on the Roe’s averaged 15 variables from the left and right nodes. Quantities with subscript 
M utilize Roe averaged values unless otherwise noted. Ideally, the difference terms used to compute the 
upwind dissipation would be in the direction of the unit normal (as occurs in an orthogonal, structured grid 
system or in an unstructured grid system employing equilateral triangles). An orthogonal difference in the 
general case would have to be interpolated from host or neighbor cells. The extra work of interpolation is 
avoided by accepting the difference q R — q r in the edge direction. First-order accurate flux is computed in 
Eq. 1 with dqiim = 0 and using Eq. 2 to compute dqM- Second-order accurate flux is computed in Eq. 
1 by forming a second-order dissipation (modeled on the Symmetric Total Variation Diminishing (STVD) 
scheme 16 ) defined as a difference of characteristic differences using dq^ m defined in Eq. 2. The characteristic 
differences centered at nodes L and R are constructed from the LS gradient of primitive variables evaluated 
at the nodes. These LS gradients are a function of all nodes connected by an edge - a domain of dependence 
which is somewhat in conflict with the original STVD formulation in that it uses dc^M to compute dcR. 
The scheme will revert to first-order with dq^ m = 0 (usually in the vicinity of shocks) when entries to the 
minmod function are of opposite sign. Otherwise, the minmod function returns the argument of smallest 
magnitude. 


dqM = Rm (<Tl — q r) 
dq L = R L ds-\7q L 
dq R = H R ds- Vq# 

ds = (x L - x R )i + (y L - y R )j + (z L - z R )k 
dcfoim = 20 M minmod [dq M , dq R , (dq L + dq R ) /4] 


(2) 


2. Shock Sensing Switch 

Note the use of a shock sensing switch (\>m in the definition of dq^ m in Eq. 2. It is required because 
the STVD limiter as approximated here is not sufficiently reliable; it admits catastrophic undershoots in 
temperature in the vicinity of strong shocks. The switch at face M is defined as the geometric average of the 
nodal values in Eq. 5. The transformation defined in Eq. 6 sets the minimum ( (fimin ) and maximum (< Pmax ) 
pressure ratios. The switch is keyed to the maximum pressure ratio computed across every edge emanating 
from a node as defined in Eq. 3. This definition implies that a high pressure ratio (<£m > ipmin) in any 
edge connected to a node is an indicator that the second-order formulation for conservation in the cell is 
suspect - even if it admits a non-zero value for dq^ m across a face. The switch is CO continuous between 
these minimum and maximum values. The final transformation in Eq. 7 makes the switch Cl continuous. 
Such a switch has not been required in the STVD formulation of flux on a structured grid; suggesting that 
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improvements to the STVD formulation are yet required in the unstructured context. In the present work 
( ( Pmim Pmax) = (2, 3) or (10, 20). The first set of values engages the switch earlier than the second set with 
consequence that temperature undershoot is more rigorously suppressed but spanwise variation of heating 
can be made worse. As grid adaptation is applied, the sensitivity to switch magnitude is decreased. 
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3. Entropy Fix as Function of Cell Reynolds Number 

The diagonal matrix A m in Eq. 1 contains the eigenvalues of the Jacobian matrix dfu/dq^M- The eigenvalues 
A i equal Um and Um A cm where Um is the velocity component normal to face M and cm is the speed of 
sound. Suppression of a carbuncle 20,21 in the stagnation region is achieved by limiting the minimum value 
of eigenvalues; the magnitude of the upwind dissipation is proportional to the magnitude of eigenvalues 
and non-physical solutions can arise if the dissipation is not restored. The limiting (sometimes called an 
entropy fix 22 ) is defined in Eq. 8 in which the eigenvalues are limited by a recommended user specified factor 
(0.05 < Aq < 0.5) times the Roe averaged sound speed at face M, cm- 


~ f I | Ai| > 2 A re f 1 

\ 4A r 'e/ 4” Are/ |Ai|>2A re / J (8) 

Are/ = Aq^MCm 

The factor ipM in Eq. 8 is used to scale down the entropy fix in regions where it may adversely effect 
the computation of mass, momentum, and energy transport by viscous dissipation. It is set to 1 for inviscid 
flow. A detrimental side effect of the entropy fix is to artificially increase computed levels of surface heating 
and skin friction. The effect is easily removed in the case of structured grids by reduction of the limiter in 
computational directions orthogonal to the surface (on faces parallel to the surface) in the boundary layer. 
In the case of unstructured grids, the same reduction can be realized by identifying a representative length 
(Im in Eq. 9) suitable for defining a local cell Reynolds number across face M ( ReM in Eq. 10). The 
arithmetic average of volumes, £2, on either side of face M divided by the area of face M provides a better 
metric for distance in the case of high aspect ratio cells then the absolute distance between nodes L and R. 
In essence, the desired characteristic length is the distance between cells in a direction normal to the face 
and not in the edge direction, which may be orders of magnitude larger in the case of high aspect ratio cells. 
The scale factor ipM is defined in Eq. 11 such that the entropy fix is fully engaged across faces where the 
cell Reynolds number exceeds a target value (Retarget = 500 used herein). It is exponentially reduced as cell 
Reynolds numbers of order 1 to 10 are approached. Note that cell Reynolds numbers of order 10 or less are 
required to adequately resolve boundary layers and free shear layers. 


Im = 


ReM = 


+ &R 

2 Am 
PmCmIm 


Pm 


'ipM = min 


1 , exp I 1 — 


— y 


Re 


( 9 ) 

( 10 ) 

( 11 ) 


6 of 19 


American Institute of Aeronautics and Astronautics Paper 2007-3960 



4- Supplementary Flux Limiter - Experimental 


A supplementary limiter has been tested that shows some promise for improving the quality of computed heat 
transfer. It is not yet coded as an option in FUN3D because there are still stability problems associated with 
its implementation. These mixed results are presented here to capture these ideas and perhaps stimulate 
additional development. The basic idea addresses the observation that the STVD limiters designed to 
prevent new maxima and minima from being introduced into the characteristic variable variation between 
nodes L and R do not prevent new maxima and minima from being introduced into the flux distribution. 
Consequently, Eq. 1 is rewritten as Eq. 12 and gTm is further limited using Eqs 13-15 to enforce either 
fk < fM < fp or fz, > fM > fe- The limiting factor for conservation law i simply reduces the magnitude 
of the upwind dissipation term dfi. It does not change the order of accuracy of the reconstructed inviscid 
flux. 



+ f R — d^M | 


(12) 


di M = Rjv/IAmI [dqM ~ dqum] 


(13) 


dfi = ujidfi 


(14) 


uJi = min 


1 | fi,L fi,R | 

I dfi\ +e 


where e = 10 10 is a small parameter to prevent divide by zero. 


(15) 


5. Post-Processing Diffusion Flux to Surface 

Diffusion of mass (gas-surface interactions), momentum (shear), and energy (heating) to a wall is evaluated 
directly from the appropriate conservation law applied to the median dual control volume for a node O on 
the wall as shown in Fig. 4. Inviscid and viscous flux across interior cell walls (shown in red connecting 
centroids C to midpoints of edges M) are computed in the usual manner. Inviscid flux across the green 
surface (segment Mi - Mf) centered at node O is evaluated as a function of nodal values at O and edge 
connected neighbor nodes on the surface. For a solid (non-porous) surface, only surface pressure contributes 
to a force in the momentum conservation laws. The diffusion flux across the wall is then calculated to force 
the sum of all inviscid and viscous flux across all walls to equal the volumetric source (if any) in the control 
volume. The heating rate, for example, is simply the remaining energy diffusion across the green surface 
(segment M\ - Mf) at O divided by the area of the green surface as represented in Fig. 4. This approach 
is simpler than a post processing step that backs out gradients normal to the surface using an appropriate 
stencil of interior nodes. All of the pieces of the computation are available during the setup of the residual 
and are captured in the boundary condition module for subsequent use in post-processing the solution. 



Figure 4. Representation of control volume at wall used to define surface heating rate from completion of 
energy balance. 
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IV. Hypersonic Flow over Cylinder — A Simple Test for Heating 


Accurate simulation of stagnation region heating in hypersonic flows is a key requirement for acceptance 
of any algorithm proposed for use in aerothermodynamic analyses. A structured grid solution generated with 
LAURA is used as both a benchmark and to generate initial grids for use in FUN3D. The test problem 1,2 
uses Uqo = 5000m/s, = 0.001kg/m 3 , and T ^ = 200K. Sutherland’s law for air is used to define transport 

properties in all perfect-gas cases. As noted below, this simple problem provides insight into the ability 
of a scheme to cleanly capture the bow shock, smoothly resolve the post-shock stagnation region flow and 
predict a smoothly varying heating distribution around the stagnation point. These flowfield characteristics 
are particularly sensitive to the inviscid flux reconstruction algorithm and problems that are not evident 
in well aligned, structured (hexahedral) grids are exposed in the unstructured (tetrahedral) environment. 2 
The goal of these tests is to develop a flux reconstruction algorithm that is applied identically to structured 
and unstructured grid environments with minimal discrepancies between the simulations and to test grid 
adaptation algorithms that are required to address any remaining accuracy issues. 

A. Grids 

The structured grid, adapted to the shock and boundary layer, from LAURA is tested in FUN3D. The 
structured grid is then converted from hexahedral elements to tetrahedral elements by adding diagonal edges 
consistently from minimum index to maximum index corners. A comparison of the grids in a plane of nodes 
perpendicular to the cylinder axis is presented in Fig. 5. The structured grid has 65 nodes from the body to 
the inflow boundary ahead of the bow shock and 61 nodes from the left to right outflow boundaries. Node 
placement is identical between the two grids. The placement of additional edges in the unstructured grid 
is the only difference. The strong biasing of diagonals in the grid is an intentional characteristic to expose 
algorithm deficiencies that may otherwise be averaged out in the simulation. 



Figure 5. Structured grid (LAURA) and biased, unstructured grid (FUN3D) in plane orthogonal to cylinder 
surface. 

A key element of this test problem is the addition of ten spanwise cells, shown in Fig. 6, across the 
cylinder, providing additional degrees of freedom in the simulation to allow asymmetries to develop. Earlier 
tests 2 (when the code was referred to as High Energy Flow Solver Synthesis - HEFSS) have shown that the 
single spanwise cell grids show good agreement with the structured code results in heating. (Fig. 7(a)) The 
strongly biased grids produce a streamwise asymmetry in the shear distribution and overpredict the LAURA 
results throughout even in the case of a single spanwise cell. (Fig. 7(b)) The contours of heating and shear 
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should be constant (straight lines) in the spanwise direction. Asymmetries introduced using the biased grid 
are easily observed in these tests; consequently, the simulation of symmetric results is a measure of solution 
quality. 



L. 



(a) Structured, quadrilateral surface elements 


(b) Unstructured, triangular surface elements 


Figure 6. Surface mesh on cylinder with ten rows of spanwise cells. 



(a) Heating 



Figure 7. Surface heating and shear over cylinder and standard test conditions with 5-species reacting air 
model and fully catalytic wall. 
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B. Structured Grid Results — Hexahedra 


Structured, hexahedral grids are assessed in FUN3D to confirm that spanwise symmetry is preserved. A 
perfect gas case (7 = 5/3, MW = 2 = 4.25) is tested with zero eigenvalue limiting. Previous tests 

with LAURA using a single spanwise cell confirmed that the carbuncle did not form for the special case of 
7 = 5/3 . 23 The carbuncle did not form using FUN3D as well; however, a spanwise variation developed in 
the solution as shown in Fig. 8 (a) and 9(a). A constant spanwise value of heating and shear was recovered 
when the eigenvalue limiting was engaged with Ao = 0.1 shown in Fig. 8 (b) and 9(b). The near perfect 
overplotting of symbols in Fig. 9(b) for each of 11 spanwise nodes (ten spanwise cells) at each x location 
around the cylinder indicates that the baseline algorithm is behaving as expected for the structured grid 
case. These results show that the eigenvalue limiting provides dissipation in the spanwise direction (where 
Um = 0) that serves to suppress a non-physical result. It was not necessary to apply a pressure sensing 
switch in this moderate supersonic Mach number case. 




j , - _ -1 

-1 -0.5 0 0.5 1 

X 



Figure 8. Surface heating contours on cylinder computed on structured grid with and without eigenvalue 
limiting for 7 = 5/3. 



(a) Ao = 0.0 



Figure 9. Surface heating and shear distributions for each of ten spanwise rows of cells on cylinder computed 
on structured grid with and without eigenvalue limiting for 7 = 5/3. Heating is represented by blue or cyan 
symbols. Shear is represented by red or pink symbols. 
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C. Unstructured Grid Results — Tetrahedra 

The perfect gas test case described above is repeated on the equivalent unstructured grid with identical 
location of nodes. The structured and unstructured results for heating and shear are compared in Fig. 10. 
An approximate spanwise variation of ±10% about the mean is observed in peak values of heating and shear. 
There is approximately 10% difference in the peak values of shear on the right and left sides of the cylinder 
for the unstructured grid. The spanwise mean value of the unstructured heating variation appears to be 
in good agreement with the spanwise constant, structured grid results. The spanwise mean value of the 
unstructured shear is about 20% to 25% higher than the structured grid result. These results suggest that 
the development of the velocity profile is more sensitive to the strong bias in the grid than the development 
of the enthalpy profile. It was not necessary to apply the pressure sensor switch in these moderate supersonic 
Mach number cases. 

Note that according to the criteria developed by Barth for critical Mach number as a function of 7 for 
shock instability 23 (and confirmed by the structured grid result of the previous section) the present test case 
does not produce a carbuncle. Nevertheless, the current case shows a spanwise variation in heating which 
indicates this behavior is not solely a function of the tendency of the approximate Riemann solver to admit 
a carbuncle. 




Figure 10. Surface heating and shear distributions for each of ten spanwise rows of cells on cylinder computed 
on structured grid and on biased, unstructured grid for 7 = 5/3. Heating is represented by blue symbols. Shear 
is represented by red symbols. 


The large value of 7 = 5/3 and low molecular weight MW = 2 produce a relatively large sound speed 
and moderate supersonic Mach number M ^ = 4.25. Selecting a perfect gas with 7 = 7/5 and MW = 28 
yields a lower sound speed and higher free stream Mach number M^ = 17.34. The shock capturing process 
now requires a shock sensor to prevent catastrophic undershoots in the temperature. In this case (p m i n = 2, 
( Pmax — 3, and Ao = 0.1. It is immediately evident that the spanwise variation in heating and shear shows 
significant increase in Fig. 11(a) with ±30% about the mean - an unacceptable result for aerothermodynamic 
analyses. The LAURA result in this case is 52 W/cm 2 which is in good agreement with the spanwise mean 
of the heating at the stagnation point. Fig. 11(b) shows an overlay of the Mach number contours in the 
front and back planes. It is evident that the captured shock is offset by one cell between the two planes. 
The surface heating contours in Fig. 11(c) indicate that a slight drift in the velocity - probably induced by 
grid bias - produces a higher heating toward the front (y = 0) plane. The jags in the temperature contours 
in Fig. 11(d) indicate a small overshoot in temperature in the post shock region. 

Using a metric of spanwise variation of heating, the present hypersonic result is approximately a factor 
of three worse than the previous supersonic result. Some of this difference is associated with use of the 
pressure sensor switch — some improvement has been noted if the switch sensitivity is dialed up (< fmin ± 10). 
However, a sensitivity to Mach number (through a variation of 7 and molecular weight) is still evident in 
these results. 
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(a) Heating and shear distribution for each of ten rows of 
spanwise cells. 


(b) Mach contours in front and back planes varying between 
0 and 1 in increments of 0 . 1 . 



(c) Surface heating contours. 



(d) Temperature contours in front and back plane. 


Figure 11. Simulation results for 7 = 7/5, pressure ratio switch = (2,3), and no grid adaptation. 


D. Unstructured Grid Results with Adaptation 

A grid adaptation algorithm 2 using (fM to drive adaptation is applied to the previous solution. In this 
algorithm no new nodes are introduced nor is any edge swapping enabled. The intent here is to cluster nodes 
around the captured shock to test if a cleaner shock capturing will improve the heating simulation. The grid 
across the boundary layer is left undisturbed by setting a parameter to stop grid movement if Im is less than 
the spacing at the boundary layer edge. 

Adapted grid results are shown in Fig. 12. Spanwise variation about the mean for heating is improved at 
±20%; however, this result is still considered unacceptable for aerothermodynamic analyses. Mach contours 
are almost perfectly overlayed between the front and back plane and the shock standoff appears to be equal 
on both planes when viewed on this scale. The post shock overshoot in temperature contours is less evident 
as compared to the unadapted case. Resetting the pressure sensor switch for the second-order upwind 
dissipation from ((/? m i n , (fmax ) = (2,3) in Fig. 12 to (10,20) in Fig. 13 makes no significant difference in this 
simulation. (Some improvement has been observed in other cases not reported here but the scheme is less 
robust.) 
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100 


100 



(a) Heating and shear distribution for each of ten rows of 
spanwise cells. 



(b) Mach contours in front and back planes. 


L, 



(c) Surface heating contours. (d) Temperature contours in front and back plane. 

Figure 12. Simulation results for 7 = 7/5, pressure ratio switch = (2,3), and grid adaptation. 
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(a) Heating and shear distribution for each of ten rows of 
spanwise cells. 



(b) Mach contours in front and back planes. 



(c) Surface heating 

Figure 13. Simulation 


contours. (d) Temperature contours in front and back plane. 

results for 7 = 7/5, pressure ratio switch = (10,20), and grid adaptation. 
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Detailed views contrasting the unadapted grid (Fig. 14) with the adapted grid (Fig. 15) clearly show 
that the captured shock is offset between the front and back planes but the adaptation compresses the grid 
around the captured shock and drives the position of the captured shock in the bounding planes more closely 
together. 



(a) Adapted mesh in front plane with pressure contours indi- (b) Adapted mesh in back plane with pressure contours indi- 
cating location of bow shock where it crosses the stagnation eating location of bow shock where it crosses the stagnation 

streamline. streamline. 



(c) Overlay of front and back plane grids with pressure con- (d) Overlay of front and back plane grids with focus on near 

tours. wall, stagnation region. 

Figure 14. Original grid for 7 = 7/5 and pressure ratio switch = (2,3). 
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(a) Adapted mesh in front plane with pressure contours indi- (b) Adapted mesh in back plane with pressure contours indi- 
cating location of bow shock where it crosses the stagnation eating location of bow shock where it crosses the stagnation 

streamline. streamline. 




(c) Overlay of front and back plane grids with pressure con- (d) Overlay of front and back plane grids with focus on near 

tours. wall, stagnation region. 

Figure 15. Adapted grid for 7 = 7/5 and pressure ratio switch = (10,20). 
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E. Application of a Supplementary Flux Limiter 

The supplementary limiter of Eqs. 12 - 15 in which the reconstructed flux value is limited to lie between the 
respective left and right states was applied to the converged state of the previous, adapted solution. Results 
of this application are presented in Fig. 16. Spanwise variation in heating has improved from ±20% to ±10%. 
A slight jag in the heating distribution at the stagnation point suggests that the smoothing effect of the 
entropy fix has been partially undone. An irregularity just downstream of the sonic line is also evident. The 
limiter, as currently implemented, is unsuitable for use in the flow solver. However, the improved heating 
suggests that attention be addressed not only to limiting of the characteristic variable differences but also 
to the flux differences. 




(a) Heating and shear distribution for each of ten rows of (b) Mach contours in front and back planes, 

spanwise cells. 



(c) Surface heating contours. 


(d) Temperature contours in front and back plane. 


Figure 16. Simulation results for 7 = 7/5, pressure ratio switch = (10,20), grid adaptation, and supplemental 
flux limiter. 


V. Concluding Remarks 

Computation of stagnation region heating in hypersonic flow on tetrahedral grids is the central focus of 
this paper. An argument is made that a robust scheme for aerothermodynamic analyses must either: (1) be 
free from any requirement of special topological grid constraints in resolving all flow structures; or (2) be 
able to automatically adapt specially required grid topologies wherever they are needed in the flow field. It 
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is noted that three-dimensional simulations using high aspect ratio tetrahedra to resolve the boundary layer 
and the bow shock produce poor heating. This problem is more severe as the Mach number is increased. 
Consequently, standard practice has evolved to use special topological elements (prisms) aligned with the 
shock and boundary layer to overcome the greatest source of difficulties in hypersonic heating simulations 
and accept the errors (or the need for a finer grid) that are necessarily introduced in such practice. 

The position advocated herein is that unstructured grids provide the greatest flexibility to adapt to 
evolving flow structures (viscous and inviscid) and complex, deforming bodies in a hypersonic flow simulation 
without a requirement for significant user intervention. It provides the greatest opportunity to create a robust 
aerothermodynamic simulation. A discussion of the simulation needs for hypersonic flow over an inflated 
aerobrake in the provide context for this position. The problem with the simulation of heating on tetrahedral 
grids is currently one of the biggest obstacles to realization of such a simulation capability. 

The goal of this paper is to ultimately define an algorithm which produces consistently accurate results 
regardless of the topology of the cells. Several algorithm modifications to handle strong shocks in the context 
of a quasi-one-dimensional reconstruction on tetrahedral grids are discussed. Evaluation of the diffusion of 
mass, momentum, and energy to the surface with a conservative formulation is defined. Nevertheless, a 
robust solution to the problem has not yet been identified. The test cases indicate: 

1. high aspect ratio, tetrahedral cells produce good heating when a single spanwise cell is utilized and 
spanwise velocity drift is suppressed; 

2. even when the entropy fix can be turned off for cases in which carbuncle formation is not an issue a 
deficiency in heating predictions on the ten spanwise cell test case is evident; 

3. spanwise variation of heating can be admitted even for structured grid simulations if the entropy fix is 
suppressed; 

4. grid adaptation at the shock with no grid movement in the boundary layer makes significant improve- 
ment to heating quality; 

5. a supplemental limiting process on the flux difference in addition to the characteristic variable differ- 
ence reduces spanwise variation but introduces some low level instabilities behind the captured shock 
downstream of the sonic line; and 

6. the ten spanwise cell cylinder test case with biased grid is a very sensitive indicator of algorithm 
robustness in that the constant biasing will not average out any deficiencies in the shock capturing 
algorithm. 
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